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We investigate the Josephson effect between two coupled superconductors, coupled by the tun- 
neling of pairs of electrons, in the regime that their energy level spacing is comparable to the bulk 
superconducting gap, but neglecting any charging effects. In this regime, BCS theory is not valid, 
and the notion of a superconducting order parameter with a well-defined phase is inapplicable. 
Using the density matrix renormalization group, we calculate the ground state of the two coupled 
superconductors and extract the Josephson energy. The Josephson energy is found to display a 
reentrant behavior (decrease followed by increase) as a function of increasing level spacing. For 
weak Josephson coupling, a tight-binding approximation is introduced, which illustrates the physi- 
cal mechanism underlying this reentrance in a transparent way. The DMRG method is also applied 
to two strongly coupled superconductors and allows a detailed examination of the limits of validity 
of the tight-binding model. 
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I. INTRODUCTION 

The Josephson effect, i.e. the flow of a zero-voltage 
current between two weakly coupled superconductors, 
with a sign and amplitude that depends on the differ- 
ence of the phases of their respective order parameter, 
can be regarded as one of the most striking illustrations 
of phase coherent behaviour in a macroscopic system and 
as one of the hallmarks of superconductivity. Although 
the Josephson effect is in general well understood, there 
is still a regime in which it has not yet been studied in 
detail: superconductors that are so small that the dis- 
crete nature of their energy levels becomes important. 
In this regime, the theory of Bardeen, Cooper and Schri- 
effer (BCS), which the quantitative understanding of the 
Josephson effect has been based on, is not applicable. 

This is because one of the underlying assumptions in 
standard BCS theory is the presence of a (quasi-) con- 
tinuous energy band. As Anderson first pointed out0, 
BCS theory is not consistent anymore once the super- 
conductor is so small that the mean level spacing d is of 
the order of the superconducting gap Abcs : Accord- 
ing to BCS-theory, the dominant contribution to pairing 
correlations comes from levels within a range of order 
Abcs around the Fermi surface, but there are no levels 
left within this range when d > Abcs ■ 

When it became possible to reach this regime ex- 
perimentally by doing transport measurements on su- 
perconducti ng gr ains with a diameter of only a few 
nanometers [j, [Ty], interest was spurred in a description of 
the pair-correlated state that is also valid for d > Abcs ■ 
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It turned out that the BCS interaction in this regime 
had already been extensively studied in the context of 
nuclear physics, where an exact solution of the reduced 
BCS Hamiltonian with discrete energy levels had been 
found in 1964 by Richardson[l5j. 

Using this solution, it was possible to explore in de- 
tail the breakdown of BCS theory as d increases. Sev- 
eral surprising insights were gained, one of which be- 
ing that BCS theory already becomes unreliable when 
d > A| cg /wDobyo, in other words, long before the An- 
derson criterion is met|l6|. The underlying reason is that 
in this regime, BCS theory underestimates the contribu- 
tion of the so-called "far or distant levels", i.e. energy 
levels farther away than Abcs from the Fermi surface. If 
the contribution of these levels is properly accounted for, 
remnants of superconductivity turn out to persist even 
for d > Abcs- Indeed, the recent experiments on small 
superconducting grains 0,^1 indirectly confirmed these 
results, in the sense that even for level spacings as large 
asd~ Abcs > they observed an even-odd effect as a clear 
indication of remaining superconducting correlations. 

Another issue that arises for small superconductors is 
that the superconducting phase (j) is not well-defined: 
When the mean number of electron pairs (TV) is so small 
that fluctuations around (N) in the grand canonical en- 
semble are not negligible anymore, N has to be treated 
as fixed. As a consequence, due to the uncertainty rela- 
tion [N, <f>] = i, the notion of an order parameter with a 
well-defined phase loses its meaning. 

Therefore, a very natural question arises: What is the 
fate of the Josephson effect between two small super- 
conducting grains, in a regime where BCS theory breaks 
down, and where the notion of a superconducting phase 
variable is no longer valid? 

In this paper, we examine this question in detail by 
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studying two pair-correlated grains, coupled by a tun- 
neling term that allows pairs of electrons to tunnel be- 
tween the grains. We study this system using the density- 
matrix renormalization group (DMRG), a powerful nu- 
merical approach applicable to strongly correlated sys- 
tems, which has already proven to be useful for calculat- 
ing the properties of single superconducting grains[5ll7|. 
We here use it to calculate the ground state of two cou- 
pled grains and to extract the Josephson energy. For 
weak Josephson coupling, we also perform a tight-binding 
approximation and compare its results to those of the 
DMRG calculation for two coupled grains. 

We identify two competing effects due to the discrete- 
ness of the energy levels: Somewhat surprisingly, the 
Josephson energy is found to be enhanced for large level 
spacing due to the contribution of a single energy level. 
At intermediate level spacing, a kinetic energy term dom- 
inates, which suppresses the Josephson energy. The com- 
petition of these effects leads to a surprising reentrant 
behavior (decrease followed by increase) of the Joseph- 
son energy as a function of increasing level spacing. In 
the limit of vanishing level spacing, the BCS result is 
recovered. 

At this point, we should mention an important restric- 
tion on our analysis: In the regime of small superconduc- 
tors that we are interested in, the charging energy for an 
electron pair to tunnel between the two superconductors 
can become huge, easily of the order of a few hundred 
Kelvin in the experiments of [19|. As will be explained 
in some detail in subsection III Dl the dominant effect 
of the charging energy is to suppress tunneling events 
altogether and thereby to destroy the Josephson effect. 
However, the interest of the present paper is to study 
the effects due to the discrete spacing of the energy lev- 
els rather than that of charging effects, which have been 
thoroughly examined already j'l lii Il2|. Therefore, we 
set the charging energy to zero in this paper. 

To experimentally realize the no-charging-energy 
model studied here, one needs systems for which the 
mean level spacing is larger than the charging energy. 
In principle, it is possible to reduce the charging energy 
of isolated grains, e.g. by using a pancake-shaped grain 
geometry in order to increase the inter-grain capacitance 
area, or by embedding the grains in a strong dielectric 
medium. - A more radical and at this point purely 
speculative way of studying Josephson physics in the ab- 
sence of charging effects would be to use uncharged par- 
ticles instead of electrons, e.g. a degenerate Fermi gas 
of charge-neutral cold atoms in a double-well trapping 
potential. Although a "superconducting" phase for cold 
neutral fermionic atoms has not yet been observed, there 
are predictions that this should be possible 8] . Once this 
has been achieved, a natural next step would be to study 
the Josephson effect in this system, for which the charg- 
ing energy would indeed be zero. 

The outline of the paper is as follows: In Section II, we 
review the theory of the Josephson effect in a way that 
is also applicable for small superconductors, for which 



standard BCS theory is not applicable, and we give a 
definition of the Josephson energy independently of a su- 
perconducting phase variable. Section III contains a brief 
introduction to the DMRG method and its application to 
the system of two coupled superconductors. Finally, in 
section IV we present and discuss the results of our cal- 
culation. 



II. JOSEPHSON EFFECT FOR WEAKLY 
COUPLED SUPERCONDUCTORS: THEORY 

In this section, we review some standard results of the 
theory of the Josephson effect and explain their relation 
to the pair-tunneling models to be used below. Our dis- 
cussion of the Josephson effect is restricted to weak cou- 
pling between superconductors, such that perturbation 
theory in the coupling can be applied. We are careful, 
however, to formulate the Josephson effect in such a way 
that a generalization beyond perturbation theory is pos- 
sible; this is done in the last subsection III Fl 

The physical assumption underlying perturbation the- 
ory is that the tunnel coupling between the supercon- 
ductors is so weak that it is energetically not favorable 
to create excited states with broken electron pairs in the 
individual grains. Therefore, the low-energy states of the 
coupled system will not contain any of these excitations, 
which will only be present as virtual states in pertur- 
bation theory. In more quantitative terms, the weak- 
coupling condition is Ej <C A sp , where A sp is the lowest 
energy of a pair-breaking excitation, and the Josephson 
tunneling matrix element Ej is defined in Eq. I|I5|) below. 

We are also careful to formulate our discussion of the 
Josephson effect independently of the notion of a super- 
conducting phase variable, such that it remains valid in 
the regime of small superconductors. The material in 
this section is mostly not new and has been discussed 
in one way or the other previously 0, 0] , but we feel it 
is worth presenting it in a way that makes the ensuing 
application to small grains evident. 



A. Josephson effect as a phase dependent 
derealization energy 

In the grand canonical ensemble, the phase of a super- 
conductor cj) can be defined via the action of the pair 
annihilation operator [2lJ bi = CjfCjj_, and the state \(f>) is 
said to have a phase <f> if 

<0|6i|0> ~ e**, (1) 

4> being independent of the state i (this is the case for the 
ground state of a superconductor) . A familiar example is 
the well-known BCS ansatz wave function | (ft) — JT^ (ui + 

Vie l ^bl)\0), where it, and Vi are real. Eq. |JTJ implies that 
a state with definite phase 4> must be a superposition of 
many states |JV), each of which has a fixed number N 
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of electron pairs: 



E 

Af>0 



C N e lN *\N) 



(2) 



subject to the condition that ( JV | bf JV + 1) is real, and 
with real coefficients Cm- 

In the canonical ensemble, however, where the num- 
ber of electron pairs N is fixed, the expectation value 
Q vanishes, and the notion of a superconducting phase 
4> is obviously not valid. Nevertheless, the concept of a 
phase difference tp between two coupled superconductors 
("left" and "right", say) is still applicable, because the 
number of electron pairs on each individual superconduc- 
tor need not be definite as long as the total number on 
both superconductors is fixed. In analogy to Eq. tp 
can, then, be defined as 



{tp\b r b\\tp) 



(3) 



Here, the operators bi and b r refer to energy levels I, r 
of the left and right superconductors, respectively. As in 
Eq. 0), one has to assume that the phase in Eq. @ is 
independent of the levels I and r for tp to be well-defined. 

An example of a state with definite phase difference tp 
is, in analogy to Eq. J2J), 



N/2 

E 

=-N/2 



(4) 



with real coefficients C„. Here, the states \ v) denotes ar- 
bitrary states with N/2 — v pairs on the left and N/2 + v 
pairs on the right superconductor, subject to the condi- 



tion that (v\bj) 



r u l 



1) is real. 



We will be only interested in situations for which \v) 
has the form 



\v) = \N/2-v) L ®\N/2 + u) R , 



(5) 



where jro)/,^ are the superconducting ground states of 
the isolated L- ("left") or R- ("right") superconductors, 
each containing a definite number of pairs, n. These 
states can always be chosen to satisfy the above reality 
condition. 

As was pointed out by Josephson, the presence of a 
phase difference tp as in Eq. J3J has observable conse- 
quences when two superconductors are coupled: In par- 
ticular, for weak coupling the coherent tunneling of pairs 
induces a zero- voltage current, 

/ = I j sin tp, (6) 

that explicitly depends on tp. As is well known |H Il8|. 
the Josephson current can, via the relation 



/= [2e/h)dE/d{tp), 



(7) 



also be interpreted as a dependence of the total energy 
E on the phase difference tp. For example, we expect 
Eq. 10 to follow from the energy-phase relation 



E(tp) = const — Ej cos tp, Ej — (h/2e)Ij 



(8) 



which we will derive explicitly in III Dl in the limit d — > 0. 
A more general definition of Ej, consistent with Eq. JSJl, 
will be given in subsection III Fl where we associate Ej 
with the energy gain in the ground state (i.e. tp — 0) due 
to the coherent tunneling of electron pairs. 

The Josephson energy Ej sets the energy scale relevant 
for the Josephson effect: It is a derealization energy 
that characterizes the coupling of two materials, their 
tendency to have the same phase and the maximum su- 
percurrent Ij = (2e/h)Ej that can flow between them. 



B. Pair tunneling Hamiltonian 

Only processes that depend on the relative phase tp 
are relevant for the Josephson effect, as is illustrated in 
Eq. 0. Because of Eq. (J3J), such processes require the co- 
herent tunneling of electron pairs; therefore, they have to 
be treated at least in second order in the tunneling of sin- 
gle electrons. The main goal of this subsection will be to 
derive an effective pair-tunneling Hamiltonian, Eq. 
below, that arises at this order. 

Consider two superconductors L and R (left and right), 
each having equally spaced energy levels with level spac- 
ing d, and each with a reduced BCS interaction with 
(dimensionless) coupling constant A: 



Hi 



E 



IV 



"l, ( T l -V i ' /' ! ' 



Hji similarly, 



(9) 

where ei — Id is the bare energy of level I, a —\, | is the 
spin, and the sums are over all energy levels closer to the 
Fermi surface than the Debye energy WDebye- 

Let L and R be coupled by single electron tunneling 
with constant tunneling matrix element t, 



h.c. 



(10) 



The coupling (|10|l lowers the total energy by generating 
states such as (@J , that superimpose different numbers of 
electrons on each superconductor. For simplicity, we 
assume the sum in Eq. (I1U|I to be cut off at WDebyo in all 
numerical calculations below. 

To second order in H le , the tunneling processes can be 
described by the effective tunneling Hamiltonian 



h 2 = -j: 



H\ e rlav) {rlav \H\ e 



rlov 



(11) 



acting on the space spanned by the states | v) , defined in 
Eq. J^J. The sum in Eq. (| 1 1 1) runs over all possible in- 
termediate states \rlau) that can be reached by removing 
a single (ra)-electron from state \N/2 — v)r and adding 
a single (Zcr)-electron to state | N/2 + v) l. E t i v is the cor- 
responding excitation energy relative to the energy of the 
state We assume, however, that for given r,l,a,v, 
all states except the one with the lowest energy give a 
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negligibly small contribution to the sum, because of the 
following argument: In the BCS limit, which is valid for 
large enough values of A, all excited states are described 
by the quasiparticle operators |l8| 



where 



E°j/2 = -(v\Hj\v+1), 



(15) 



P^c 



where P^ is an operator that creates an additional pair. 
In this limit, it is easy to see that only the lowest energy 
state \rlav) = r Yr e ) rT i'~f(f l \r^ a -\ r \ 1 ') gives a contribution to 
Eq. I|l 1|) . whereas all other intermediate states have a 
vanishing overlap with Hi e \v). This is also the case for 

A = 0, where 7(e) CT ;7(fo)(- ff ) r = c^c CTr . For intermediate 
values of A, no simple argument can be made; we expect, 
however, that still the state with the lowest energy will 
give the dominant contribution. 

The energy E r i v is given by the collective excitation 
energies E r + Ei, arising from the fact that levels r and 
I are singly occupied. In general, it will also include a v 
dependent contribution from charging energy due to the 
electron tunneling, see 01 > but we chose to consider only 
situations in which these can be neglected. 

In Eq. two kinds of tunneling terms are present: 
on the one hand, terms proportional to b\b r or to b\b l 
that describe coherent pair tunneling, on the other hand, 
single electron terms proportional to c r<J c\. a c\_ lJ c l _ (T that 
describe the tunneling of a single electron from / to r 
and back. When the former terms are applied to a state 
\tp), defined in Eq. |4"|. they produce a phase dependent 
energy shift. In contrast, the latter terms only lead to 
a phase- independent energy shift, which is irrelevant for 
the Josephson effect. For this reason, the single elec- 
tron terms can be omitted from the Hamiltonian , as 
long as only phase dependent processes are of interest 0. 
Then, one finally arrives at the pair tunneling Hamilto- 
nian 



u t Pc\ a T ViC. 



Hj 



7 d 2 



E r + El 



(b\bi + h.c), 



(13) 



with 7 = t 2 . We shall use for the excitation energies their 
BCS values, E r j = JA% CS + e£,. 



C. Tight-binding model 

In the space spanned by all states without any pair- 
breaking excitations, i.e. all states of the form \ u) defined 
in Eq. (|5J), the Hamiltonian H = Hl + Hp + Hj looks 
like a tight-binding Hamiltonian: 



/ E(i 



H = 



-E" } /2 



~E°j/2 E(u) 
-E°j/2 

: 



-E°j/2 



\ 



(14) 



E{v) = (v\(H L +H 



(16) 



Much fewer electrons will tunnel than are present on 
each superconductor, v <C N. Therefore, the off-diagonal 
elements Ej can be taken to be independent of v. The 
diagonal elements are given by 



E(u) = const + 2d(v - is ) 2 



(17) 



which has the form of an effective charging energy term. 
This is because changing v by one is equivalent to shift- 
ing the relative chemical potential between the grains by 
the amount 2d, except for the outmost energy level (i.e. 
the level closest to the cutoff at ojDebye)> which can be 
neglected if v -C N. 



D. Discussion of the tight-binding model 

In this subsection, we first discuss the above tight bind- 
ing model HI 4(1 in the limit d — > and check that it is 
consistent with the well-known result of Ambegaokar and 
Baratoff pj. Then, we draw attention to what changes 
will occur as the superconductors become smaller. 

For d — > 0, the diagonal elements JHJ of Hj [Eq. ltT4l ] 
are independent of v. Also, BCS theory is valid, so the 
off-diagonal elements 1 (15( 1 are given by 



E°j 



E 



EiE r (Ei + E r ) 



2„2 



(18) 



In the left equality of Eq. 1(18(1 . the BCS expression 
for the matrix elements (y \v + 1) = viuiv r u r = 
A. 2 iCS /(2EiE r ) has been used. For the right equality, the 
sum has been replaced by an integral, = I^°oc ^d^ 2 • 
No harm is done extending the integral range beyond 
WDebyc to infinity, because it is naturally cut off at the 
scale Abcs anvwa y, assumed to be much smaller than 

^Debyc ■ 

The energy eigenstates of 1(14(1 then have the form of 
Eq. (0J with constant coefficients C v . As anticipated in 
Eq. JSJ), they correspond to an energy E(<p) = const — 
E°j cos ip, and therefore we can identify 



E, 



E°j = At 2 ir 2 



irhA 



BCS 



4e 2 R N 



(19) 



The last equality expresses Ej in terms of the normal- 
state conductance i?^ 1 = (Ane 2 /fr)t 2 , and agrees with 
the well-known Ambegaokar-Baratoff formula[l| at zero 
temperature. 

Now we turn to the question what happens when the 
superconductors enter the regime d > Ag CS /a;Dcbyc, in 
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which the BCS Ansatz wavefunction becomes inappro- 
priate ^(|. The transition to this regime is straightfor- 
ward now, because the tight-binding model itself remains 
valid: The diagonal and off-diagonal matrix elements 
E(v) and Ej of the tight-binding Hamiltonian f!4jl. de- 
fined in Eq. (|17(l and Ijl5(l , will no longer be given by the 
BCS expression, but will have to be evaluated using the 
exact ground state wave function: The effect of the dis- 
crete level spacing on the diagonal elements E(v), given 
by Eq. (|17l) . will be to lift the degeneracy among them, 
thereby suppressing pair tunneling. 

The off-diagonal elements Ej will change with respect 
to their BCS value l|18[) due to several effects when the 
superconductors become small: (i) The excitation ener- 
gies E r \ v in Eq. (|13|l will include a term from the charging 
energy in the intermediate state, as studied in |l2|. We 
neglect this effect, because we have chosen not to study 
charging effects at all. Furthermore, the change in su- 
perconducting correlations due to the finite level spacing 
will affect both (ii) the excitation energies E r and Ei 
in Eq. (|13|l (which we, however, replace with their BCS 
value) and (iii) the matrix elements that en- 

ter Ej. Finally, (iv) the shift of the Fermi level between 
the states \v) and \v + 1) will also change the matrix 
elements Ej, as is explained in section llV Al below. 

As it turns out (see section ITV Al below) . E°j increases 
with increasing level spacing d, mainly due to the effect 
(iv). Once Ej becomes comparable to A sp , the lowest 
pair breaking ("single-particle") excitation energy, the 
superconductors can no longer be considered as weakly 
coupled, and the tight-binding model itself loses its va- 
lidity. 



E. The effect of charging energy 

As mentioned in the introduction, the Coulomb charg- 
ing energy plays an important role in small superconduc- 
tors. Although we shall neglect it in the remainder of this 
paper, here we present a brief qualitative discussion of its 
main effects. The charging energy Ec = (2e) 2 /C, C be- 
ing the inter-grain capacitance, is the energy cost for tun- 
neling an electron pair from one grain to the other. It in- 
troduces an additional term in l|14l) . E(v) ~ Ec{v — vq) 2 . 
Ec can become huge in the small grain limit and essen- 
tially destroys the Josephson effect, since it suppresses 
pair tunneling. 

Even if a gate is used to make two states | v) and | v + 1 ) 
degenerate by a suitable choice of the gate voltage (i.e. 
vq = 1/2 plus an integer), such that at least one pair 
can still tunnel between the grains at no energy cost, the 
charging energy might nevertheless destroy the Joseph- 
son effect altogether: It may cause one electron pair to 
break into two unpaired electrons, one on each grain, if 
the associated lowering of the charging energy exceeds 
the energy necessary to form a pair-breaking excitation. 

An order-of-magnitude estimate shows that this ac- 
tually happens in the regime that the level spacing is 



important, if no measures are taken to reduce the charg- 
ing energy: (i) As explained above, the charging energy 
must be smaller than the lowest energy of a pair break- 
ing excitation, Ec < A sp , such that no pair breaking 
excitations occur, (ii) Abcs < \J ^Dobyc d must be sat- 
isfied if the grains are to be small enough so that devia- 
tions from BCS become important (the 'weak' criterion 
in ^3)- ( m ) For the present purpose of constructing 
an order-of-magnitude estimate, we take A sp ~ Abcs, 
although these two e nerg y scales may not be identical in 
the small-grain limit |lfi| . (They differ, for example, by 
a factor of up to two for the parameter range shown in 
Fig. I of 0.) 

Putting (i) to (iii) together, the inequality 

E c < V^Dcbycd, (20) 

independent of A, has to be satisfied. 

Let us now explore what this implies for real Alu- 
minium grains: If the inter-grain capacitance is modelled 
by an Aluminium oxide (e « 8) layer of thickness 15 
Aand area nr 2 , then Ec ~ 0.8 eV(r/nm) -2 . (A smaller 
thickness D in principle linearly decreases the charging 
energy, but at the same time, the inter-grain coupling 
t is exponentially increased [lflj. t 2 oc exp[— D/(0. 54A)]. 
Since at a thickness of less than ~ 15 A, the grains 
are so strongly coupled that thay can no longer be con- 
sidered as distinct, this distance seems to be a realis- 
tic order-of-magnitude lower bound for D.) Using the 
Debye energy WDebyc = 35 meV for Aluminum, we ob- 
tain y'wDebye • d = 0.054 eV(r/nm)- 3 / 2 , and Eq. ^ 
implies r > 250 nm. At such a large size, Aluminum 
is well in the BCS regime. According to criterion (ii) 
above, deviations from the BCS approach for a grain of 
that size would be observable only for a material with 
Abcs < 10 -5 eV, an order of magnitude less than Al. 

As an example, we consider the experiments of Naka- 
mura et al , which use a superconducting island with 
A pa 230 /ieV and Ec ~ 117 fieV. These islands are evi- 
dently so small that they are quite close (up to a factor 
of 2) to the regime where the charging energy would be- 
gin to suppress pair tunneling by favoring single-particle 
excitations. Nevertheless, their islands are still large 
enough to be well described by BCS theory. 

However, as mentioned in the introduction, the interest 
of this paper is to study the effects due to the discrete 
spacing of the energy levels, since the charging effects 
have already been discussed previously 0, |tj [2] • There- 
fore, we henceforth set the charging energy to zero. 



F. Generalization to strong coupling 

In the weak coupling limit, we have defined the Joseph- 
son energy via the part of the energy JSJ) that depends 
on if. However, Eq. (jHJ is only valid for weak coupling 
(i.e. in second order in the single electron tunneling). We 
may equivalently define the Josephson energy as the max- 
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imally possible energy lowering due to coherent pair tun- 
neling, i.e. when single electron terms are neglected as in 
the derivation of l|13|l : 



Ei 



E, 



coupled 



E, 



uncoupled- 



(21) 



This definition agrees with the usual one (JHJ in the weak- 
coupling regime, because the maximally possible energy 
lowering occurs at phase difference ip = 0. Eq. H21fl allows 
an extrapolation to strong coupling as well, and therefore 
we will use it henceforth. 

Unfortunately, the pair tunneling Hamiltonian <|13|) . 
being only derived in second order perturbation theory, 
loses its validity for strong coupling; in general, one would 
have to use the single-electron tunneling Hamiltonian 
(|10|l in that case. For simplicity, however, we choose for 
our strong coupling analysis a somewhat different cou- 
pling term that only includes pair tunneling, 
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A 



BCS 



E 

rl 



h.c. 



(22) 



and that differs from the pair tunneling Hamiltonian l|13|) 
in that the intermediate energy E r + E\ has been re- 
placed by the constant Abcs- Therefore, the Hamilto- 
nian (|22|l and i|13|) are not equivalent. It is neverthe- 
less interesting to study the Hamiltonian H22JI for sev- 
eral reasons: Firstly, it captures the essential physics of 
the Josephson effect in a simple way: two superconduc- 
tors coupled by a tunneling barrier that allows for pair 
tunneling. After all, it is the pair tunneling and not 
the single electron tunneling that is at the heart of the 
Josephson effect. Secondly, for 7d/ABcs = \ the to- 
tal Hamiltonian looks just like one single superconduc- 
tor, thus (|22|l is able to describe the transition to the 
strong-coupling regime where two superconductors effec- 
tively become one. Thirdly, it is amenable to a rather 
straightforward treatment by the DMRG approach (in 
contrast, H'j of Eq. i|13|) would require much more nu- 
merical effort), which has the very significant advantage 
of yielding direct access to the regime of strong coupling 
between the two superconductors. 

At weak coupling, a tight-binding analysis for (12211 sim- 
ilar to the one that led to Eq. Ijl8(l can be performed. In 
the large grain limit, one finds the Josephson energy to 
be 



E°j = 



2 7 A 



B( 'S 



A 2 



(23) 



independent of d. In other words, the Hamiltonian 1221) 
has a well-defined continuum limit when 7 is held contant 
as d — > 0, as it should. 



III. DMRG APPROACH 

In the context of nuclear physics, Richardson found an 
exact solution .15] of the Hamiltonian © for a single 



superconductor, that allows in principle to calculate all 
of its eigenenergies and eigenstates. Because the tight- 
binding calculation for weakly coupled superconductors, 
as outlined in III CI only needs the matrix elements H15|) 
between states of a single superconductor, Richardson's 
solution is, in principle, sufficient for that case. 

However, while the eigenenergies of (0 can be calcu- 
lated with only little numerical effort using Richardson's 
solution, the computation time needed for the eigenstates 
and for matrix elements like the ones in (|15f) scales like nl 
with the number of energy levels n in the system, mak- 
ing it effectively impossible to go beyond, say, n = 12 
levels or so (more precisely, only the number n of en- 
ergy levels between i<<Fermi - ^Dcbyc and Epcrmi + ^Dcbyc 

matters). For this reason, despite there being an exact 
solution available, it is indispensable also for the tight- 
binding model to have an alternative approach at hand 
that is approximate, but manageable. Moreover, for the 
strong coupling analysis in Sec. Ill Fl that invokes the 
pair tunneling term H22fl . Richardson's solution is not 
applicable at all, so that the use of a different approach 
becomes unavoidable. 

For these reasons, we have adopted an approach based 
on the density matrix renormalization group (DMRG), 
whose power and efficiency for dealing with pair- 
correlated nanograins has been demonstrated recently^, 
Il7j . We will use two kinds of DMRG calculations: A 
single-grain DMRG for calculating the matrix elements 
(Sec. I15|) to be used in the tight-binding model at weak 
coupling (cf. Sec. Ill Cjf) . and a two-grain DMRG for the 
case of strong coupling (cf. Sec. Ill F|l . 

In this section, we first discuss some general aspects 
of the DMRG algorithm in energy space in 1111 Al leaving 
some of the more technical issues for appendix [S] In 
Sec. IIII 51 we discuss the one-grain DMRG, and turn to 
the dicussion of the two-grain DMRG in Sec. IIII CI 



A. The DMRG method in energy space 

The DMRG in its usual implementation is a real- 
space renormalization group method, and has been very 
successful for describing one dimensional many-particle 
quantum systems, such as spin chains 01- Usually, the 
Hilbert space for such systems is too large to be diag- 
onalizcd exactly on a computer. The DMRG algorithm 
allows to keep only a reduced part of the Hilbert space 
that is small enough to be tractable even on a desktop 
computer, but still sufficient to describe one or several de- 
sired states, the so-called target states (in our case, the 
ground state will be the target state). This is achieved by 
progressively increasing the chain size, adding sites one 
at a time, while only a limited number of states is kept 
at each step, those states being selected as the most rel- 
evant ones for describing the target state(s) in a density 
matrix analysis. 

Although the DMRG is mostly limited to one dimen- 
sional systems, it can be applied to three dimensional 
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ones by using the energy axis as the one dimensional 
"system" , such that the bare energy levels play the role 
of sites on a one dimensional chain. This is not always 
useful, because the interactions between these "sites" can 
be much messier than between sites in real space, the 
latter being generally local. Luckily, as will be seen, the 
BCS interaction is, although nonlocal, simple enough for 
the DMRG algorithm to be applicable. 

The DMRG builds up the system, starting from the 
low-lying energy states around the Fermi surface, which 
are the physically most important ones, and successively 
adds levels lying further and further away from the Fermi 
energy. It should be noted that this is quite contrary to 
the way usual RG calculations are performed, where high 
energy levels are integrated out, approaching the low en- 
ergy states from above. This allows these two comple- 
mentary approaches to be simultaneously applied: As 
long as not all energy levels have yet been added to the 
system, only the ones near the Fermi surface are explic- 
itly included in the DMRG calculation. The other ones, 
which will be included only at later steps, are meanwhile 
taken into account using a renormalization of coupling 
constants (as introduced in eq. (43) of 0])- 

For this purpose, the following scheme turns out to be 
numerically very efficient for renormalizing the coupling 
constants A and 7: When the i levels closest to the Fermi 
energy are included, choose the coupling constants, say 
Xi and 7^, such that the BCS band gap Aj of the current 
system equals the final value A„, where n is the desired 
final number of levels. In the DMRG for a single grain, 

A\ 1] = id/(2sinh(l/Ai)). In the two-grain DMRG, the 

(2) 

band gap is given by A- = id/ sinh(l/(Aj + 7,-<i/A n )). 
The latter is the the solution of the BCS gap equation 
with two different interaction matrix elements — \d and 
—jid 2 /An \ as in © and (|22l) . For large couplings, this 
scheme turns out to be more efficient than a perturba- 
tive renormalization of the coupling constants. At weak 
couplings, for which perturbation theory is expected to 
work, both approaches perform equally well. 

Another drastic reduction of degrees of freedom occurs 
because in the model we study, the energy levels that are 
occupied by a single electron completely decouple from 
all the interaction terms ©, 1131) and (|22[l . Because the 
creation of a singly occupied level is associated with the 
energy A sp and therefore energetically unfavorable, there 
will be no singly occupied levels in the low-energy sector 
of a superconductor, if one assumes the total number of 
electrons to be even. Due to these considerations, we can 
omit these levels from the beginning, and consider only 
the case of empty or doubly occupied energy levels 0] . 

Although the full Hilbert space is drastically reduced 
by the DMRG algorithm, it produces excellent results. 
In the case of the two-grain DMRG, the accuracy can 
be checked by comparing the condensation energy from 
DMRG to the Richardson solution, which is available for 
two specific values of the inter-grain coupling 7 in (I22II . 
namely for ^lj id/ A — A (which effectively describes 
one single, larger superconductor) and 7 = (two inde- 
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FIG. 1: Relative error between the exact result from Richard- 
son's solution and the two-grain DMRG at BCS coupling 
A = 0.4, with n = 100 energy levels per grain. The inter- 
grain coupling in the upper plot is yd/ A = A. In the lower 
plot, 7 = 0. 

pendent superconductors). The results for the two-grain 
DMRG are shown in Fig. ^ and show the following fea- 
tures: (i) High precision at strong inter-grain coupling, 
with a relative error in the condensation energies of only 
~ 10~ 7 when m = 100 states are kept, (ii) Lower, but 
still sufficient precision for decoupled grains (7 = 0): 
~ 10~ 3 for m — 300, for n — 100 energy levels. However, 
the algorithm fails at weak coupling when the number 
of energy levels n becomes large (N > 80 — 150), see 
IIII CI In this case, a perturbative calculation fsee IIIIB"|) 
becomes necessary. 

With the one-grain DMRG, the accuracy of case (i) 
is obtained. As always in DMRG, the precision can be 
systematically improved by increasing m. 

B. One-grain DMRG for tight-binding model 

If the grains are weakly coupled, the tight-binding ap- 
proach can be applied, based on the Hamiltonian Q13p. 
Here, the microscopic model only enters via the tunnel- 
ing matrix elements Ej (|15|) of the tight-binding Hamil- 
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tonian (|14|) . Although these can in principle be calcu- 
lated exactly using Richardson's solution, in practice the 
DMRG algorithm is much better suited for that task, as 
explained above. 

Assuming v <C N and using Eq. J5|, the only matrix 
elements needed for Ej are (N/2 + 1 1 fot | iV/2) for all val- 
ues of i. We evaluate these matrix elements using the 
DMRG algorithm for one single grain, as introduced in 
Refs. HEa. This requires the simultaneous knowledge 
of two ground states with differrent pair occupation num- 
bers, |iV/2) and \N/2 + 1). These states are constructed 
in a single run, as explained in appendix^ Once these 
matrix elements have been calculated, it is straightfor- 
ward to diagonalize the tight-binding Hamiltonian 114(1 . 
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C. Two-grain DMRG 

If the DMRG is directly applied to a system of two 
grains, the regime of strong coupling can be explored, 
too. For this purpose, we use the inter-grain coupling 
term (12211 . introduced in subsection III Fl The exact 
Richardson solution cannot be applied for this system 
(except for the particular value of 7 = Ad/AecSi which 
has been used in subsection IIII Al for checking the accu- 
racy of the results). 

Although the two-grain DMRG can cover the previ- 
ously unaccessible parameter region of strong coupling, 
it turns out to fail for too weak inter-grain coupling, if 
the system is large (more than, say, 80-150 or so energy 
levels, depending on the other parameters). The reason is 
that the DMRG relies on correlations between the grains 
for being able to effectively reduce the Hilbert space, 
and these correlations vanish in the limit 7 — > 0. This 
can easily be seen in the limiting case 7 = 0, in which 
the two grains L and R are completely uncorrelated and 
can, each, be described by mi independent basis vectors 
1 1) i, |mi)i and | l)n, This implies that the 
m basis vectors retained are essentially product states of 
the form €5 \j)r, and only an accuracy corresponding 
to mi = ^fm kept basis vectors per grain is attained. If 
the inter-grain coupling 7 is increased, correlations be- 
tween grains L and R quickly develop that allow to keep 
only a few dominant ones of the product states, but for 
7 = 0, and also for very small values of 7, each of these 
states is equally important, making the DMRG highly 
inefficient. That the DMRG still works even for 7 = 
if only a few (< 80 — 150) energy levels are considered, 
is due to the fact that in this case, the necessary num- 
ber of states to be kept per grain seems to be so low 
(mi ps 15) that the ground state can still be reasonably 
well approximated. 

To summarize, the two-grain DMRG works well for 
strongly correlated systems, but produces unsatisfying 
results for the case of weak inter-grain coupling. How- 
ever, this is the regime in which perturbation theory 
can be used, as described before: thus, the two-grain 
DMRG and perturbation theory are two complementary 



FIG. 2: The Josephson energy Ej in the tight-binding ap- 
proximation, based on Eq. I|14|l . as a function of the grain 
size. Ej is defined via Eq. 12 1 1 as the additional energy gain 
due to coherent pair tunneling and is normalized to the BCS 
result i?j CS in Eq. 119H . Compared to the off-diagonal ma- 
trix element Ej (dotted line), Ej (dashed line) is reduced by 
a factor of up to 2 due to the finite-size kinetic energy term 
E(v) of Eq. JT7J- The logarith mic plot in the inset shows how 
the BCS result of Eq. 11911 is recovered as d — > 0. 



approaches; their regimes of usefulness are illustrated in 
Fig- El below. 



IV. RESULTS 

A. From small to large grains: The effect of 
discrete energy levels within the tight binding model 

Fig.|21presents our results for the Josephson energy Ej, 
defined as the derealization energy l|21l) due to pair tun- 
neling, in the tight-binding approximation, calculated 
using the coupling Hamiltonian given by Eq. I|14|) . Ej 
is plotted in units of the BCS result Ej(BCS), given by 
Eq. (|19fl . The dashed line in Fig. [21 displays the Joseph- 
son energy as a function of decreasing level spacing d, i.e. 
of increasing grain size, characterized by the number of 
discrete energy levels n between ep ± c^Debye- 

While the level spacing d is varied, the parameters 
A and 7 in l|13(l are held fixed (at the values A = 0.3, 
7 = 0.05), such that the BCS value in Eq. fTJfy of Ej 
is independent of the grain size, and a well-defined limit 
d — ► exists. 

Within the tight-binding model, we observe two com- 
peting effects, to be discussed in detail below, that influ- 
ence the Josephson energy as the level spacing d increases 
(i.e. moving toward the left side of Fig- EJ: (i) On the 
one hand, the finite-size kinetic energy term E(i>) of 
Eq. (|17l) increases, which tends to reduce the Josephson 
energy Ej; (ii) on the other hand, the off-diagonal ma- 
trix element Ej in Eq. l|14f) increases (as shown in Fig. 
121 dotted line), which tends to increase Ej. The com- 
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bination of these two tendencies leads to the reentrant 
behaviour seen in Fig. [21 particularly in the inset, with 
a remarkable increase in Ej when d becomes sufficiently 
large. 

The kinetic term (i) was discussed in section III CI We 
chose vq in Eq. I|17l) as v$ — 1/2, such that the two 
lowest-lying states are degenerate and at least one pair 
can tunnel at no energy cost between these states no 
matter how large d is. For this case, the total reduction 
of the Josephson energy due to the finite-size kinetic term 
amounts to a factor of at most 2, even for very large 
level spacing d. This is because for d — > oo, all but 
the lowest two states "freeze out", so that the tight- 
binding Hamiltonian (|14|l effectively reduces to a two- 
level system for the states \v = 0) and \v = 1), whose 
tunnel splitting is Ej/2 (hence the reduction by a factor 
of 2). Nevertheless, the reduction in Fig.[5]is seen to be 
considerable even for fairly large grains (still of order 20 
% for Ascs/d ~ 100, corresponding to n ~ 3000 levels), 
because it depends on the ratio dj Ej, where E j typically 
is a small number itself. For d <C Ej, the asymptotic 
behaviour Ej = Ej(l - v / 2d/E°) (thin dashed line in 
the inset of Fig. EJ) is found in analogy to the treatment 
of small charging energies in section 7.3 of 01 , by using 
an Ansatz wave function given by Eq. (@J) with ip = and 
C„ of Gaussian form. This Ansatz wave function turns 
out to be asymptotically correct for d « Ej |is[ . 

Next, we discuss the increase of Ej in the small-grain 
limit (ii). It is due to the fact that the matrix elements 
(N L \ b] \ N L + 1)(N R + l\bl\N R ) that contribute to E° } in 
Eq. I|15|) have a different number of electron pairs in the 
states acting from the left and on the right. This fact is 
neglected in standard BCS theory, where the total num- 
ber of pairs is assumed to be macroscopically large any- 
way. When the level spacing d becomes large, however, 
this is the main reason for the increase of Ej-. 

The increase of Ej is easy to understand for the Fermi 
state (A = 0) and in the BCS limit (A > 2/ In N, see HU). 
In the Fermi state, the matrix element (A|&i|A) is zero 
for all values of i, but (N \bi\ N+l) gives a contribution of 
1 for the one level i = in that is below the Fermi surface 
of | N + 1) and above the Fermi surface of \N). In the 
BCS case, the matrix element is given by (A|fei| N+ 1) = 

ufv^ +1 . The upper indices on u and v indicate the 
total pair occupation numbers with respect to which they 
are taken, with the effect that v N+1 has the chemical 
potential shifted upwards with respect to v N by the level 
spacing d. Thus, the product ufv^ +1 becomes larger as 
the level spacing d increases, as is illustrated in Fig. [3] 
We shall call this modification of the BCS calculation the 
"finite-d" BCS calculation. 

In Fig. |31 the finite-d BCS matrix elements (solid line) 
are also compared to the exact values obtained using 
the DMRG (filled dots). The comparison shows that 
for the levels close to the Fermi energy (i.e. the central 
level iff and the next, say, 2 levels), the finite-d BCS 
result overestimates the pairing correlations: the (quasi- 
)exact DMRG solution is seen to have a more pronounced 
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FIG. 3: Various approximations for the matrix element 
(TV 1 6^ | TV + 1} are compared. The product of the BCS co- 
herence factors u N v N (BCS, thin dashed line) is compared 
to u N v N+1 (timte-d BCS, solid line) in a grain of n = 20 
levels. Because v N+1 in the latter product is shifted to the 
right by an amount of d with respect to v , the finite-d curve 
must obviously be larger than the BCS curve for all values 
of i. Also shown (filled dots) are the exact matrix elements 
{N\bi\N + 1), calculated using the DMRG approach. The 
inset compares the BCS, the finite-d BCS and the DMRG re- 
sults for a larger grain with n = 200 levels. While the finite-d 
BCS result shows excellent agreement with the DMRG result 
in this regime, the BCS result still deviates from the exact 
result. 

Tight binding matrix element for^.= 0.3, Y= 1 
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FIG. 4: E" is evaluated in the BCS, the finite-d BCS and 
the DMRG method. Although BCS theory is not valid in the 
regime of large level spacings d, the finite-d-BCS reproduces, 
at least qualitatively, the correct behaviour seen in the DMRG 
curve. 



peak at the central level zjv, whereas the contribution 
of the neighbouring two levels is somewhat reduced, re- 
sembling, for these levels, qualitatively more closely the 
A = case discussed above. For the energy levels further 
away from the Fermi energy than that, the finite-d BCS 
calculation is seen to slightly underestimate the matrix 
elements. This is not unexpected, because BCS theory is 
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known 16] to underestimate the superconducting corela- 
tions of energy levels much farther away from the Fermi 
surface than Abcs , which for the parameters of Fig. is 
Abcs » 0.7. 

-Ej, however, being a weighted sum over all products of 
these matrix elements, is nevertheless not so far off in the 
finite-d BCS approach even for very small grains, as can 
be seen in Fig. 0] This is surprising and somewhat for- 
tituous, since the BCS theory does not self-consistently 
describe the grains in the limit that they are small. The 
reason that the finite-d BCS works so well seems to be 
that the underestimation of the matrix elements for level 
i n and and their overestimation for the other levels can- 
cel each other to a large degree. 

In conclusion, the main reason why Ej increases as 
the grains become small is very simple: the chemical po- 
tential of the grains shifts due to the finite level spac- 
ing whenever a pair tunnels from one grain to the other. 
Note that the BCS ansatz without taking this effect into 
account is not accurate near the Fermi energy even for 
fairly large grains (see the inset in Fig. |3J), for which 
the finite-<i BCS theory agrees perfectly with the DMRG 
result. 

The competition between the finite-size kinetic term 
on the one hand and the increase of Ej on the other 
leads to the reentrant behaviour of Ej as seen in Fig. 
|5] This is one of our main results. Two regimes can be 
distinguished as a function of Abcs/^: For very small 
grains (Abcs/^ < l)j superconducting correlations are 
only weakly present, but the 1-level effect outlined above 
leads to a strong enhancement of Ej and, therefore, of 
the Josephson energy Ej. Despite not being a well- 
justified approximation in this regime, the finite-c? BCS 
result nevertheless gives a surprisingly good estimate of 
the Josephson energy. On the other hand, for larger val- 
ues of Abcs/^ (> 10, say), Ej is almost constant and 
very close to the BCS value. The kinetic energy term 
in Eq. (|17fl . however, reduces the Josephson energy be- 
low the BCS value, and vanishes only rather slowly. The 
reentrant behaviour of Ej occurs at the intermediate re- 
gion 1 < Abcs/^ < 10, in which both effects are com- 
peting, and in which Ej is slowly approaching its BCS 
value from above. 



B. Limitations of the tight-binding approach 

The tight-binding approximation, which neglects pair 
breaking excitations of the individual grains, is valid only 
for small couplings, such that Ej lies well below the low- 
est excitation energy A sp . In Fig. however, Ej is seen 
to grow strongly with increasing level spacing d. Thus, 
for sufficiently large d, the tight-binding approach in- 
variably becomes unreliable, and a different method is 
needed. In order to complement the tight-binding ap- 
proach and to check its quality, we have thus used the 
two-grain DMRG solution that does not rely on the inter- 
grain coupling being weak. 
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FIG. 5: A rough sketch of the regimes of validity for the 
DMRG and the tight-binding approach in parameter space 
(inter-grain coupling 7 vs. the number n of energy levels). 
There is only a small overlap (shaded region) at small 7 and 
small n, in which both approaches simultaneously work well. 



The DMRG, however, has its own limitations, as was 
explained in subsection IIII CI Firstly, it requires a pair 
tunneling Hamiltonian H22[) that describes a somewhat 
different model. This implies, of course, that it has to be 
compared to a tight-binding model using the same pair 
tunneling Hamiltonian as well. Secondly, the two-grain 
DMRG can break down at small couplings if the number 
of energy levels is large, for precisely the same reason that 
the tight-binding model works well: The correlations be- 
tween the two grains, which the DMRG relies on, become 
very weak. 

The regimes of validity of the two complementary 
approaches are schematically depicted in Fig. [5] The 
tight-binding method only works well at small coupling, 
A s . p . <C Ej (region left of dashed line), whereas the 
DMRG works well only at large coupling (region right 
of solid line). A simple (analytical) condition for the va- 
lidity of the DMRG approach cannot be given, which is 
why the axes in Fig.JSJare drawn without units. However, 
the quality of the DMRG approach is found to depend 
sensitively on the number of energy levels n. In particu- 
lar, the DMRG turns out to be reasonably accurate for 
all values of 7 down to 0, as long as n < 80 — 150 (depend- 
ing on other parameters), as is motivated in subsection 
HIEO and seen in Fig. Eland Fig.Q 

In Fig. [5] and Fig. [7| the tight-binding approximation 
for the Josephson energy is compared to the two-grain 
DMRG as a function of the grain size, for two differ- 
ent values of the inter-grain coupling 7 (corresponding 
to moving along vertical lines in Fig.[SJ). The Josephson 
energies are again plotted in units of their BCS values, 
now given by Eq. (jSSJ- In Fig.© both methods are seen 
to agree for small numbers of energy levels, n < 80—100. 
For larger values of n, the two-grain DMRG breaks down, 
for the reasons outlined in lHTCl The DMRG method it- 
self signals its own breakdown: Convergence as function 
of the kept DMRG states m is no longer achieved, as can 
already be seen when the two curves shown in Fig. [5] 
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FIG. 6: The Josephson energy is calculated in the tight- 
binding and in the DMRG approach. In agreement with Fig. 
|SJ both curves agree if the number of energy levels n is small, 
but the DMRG fails for n > 80 - 100. 



which correspond to m — 330 and m = 360, are com- 
pared. 

Since both the two-grain DMRG and the tight-binding 
approach are ultimately variational methods, the one 
that produces the higher value of Ej (i.e. the lower total 
condensation energy) must be the better approximation. 
Also in this respect, the DMRG method is seen to be 
failing for n > 80 — 100 in Fig. EI m agreement with Fig. 

El 

Fig. [7| shows the result of a similar calculation as Fig. 
EJ at a higher value of the inter-grain coupling 7 = 0.01. 
Now, the two-grain DMRG is seen to be valid up to some- 
what larger values of n. For small n, n < 100, the DMRG 
now produces a higher value of Ej , indicating that in this 
regime, it produces a better result than the tight-binding 
method, as anticipated in Fig. [SJ 

The results from Fig. |S] and Fig. [7| are similar to the 
ones in Fig. [21 where only a tight- binding calculation had 
been performed. In particular, the two-grain DMRG re- 
produces the increase in Ej for small values of n, corre- 
sponding to large level spacing d, and thereby confirms 
the reentrant behaviour observed in the tight-binding 
aproach (cf. Fig. 01. 

In Fig. |S1 the tight-binding and the two-grain DMRG 
results are plotted as a function of the inter-grain cou- 
pling 7 (corresponding to moving horizontally in Fig.|SJ. 
The plot is extended to very large values of the inter- 
grain coupling 7 in order to show the point at which the 
two-grain DMRG can be compared to the exact result at 
(c£/Abcs)7 — \ which it reproduces nicely. We empa- 
size that in the regime of large 7, some of the physical 
assumptions (e.g. the use of a tunneling Hamiltonian) of 
our calculation are not justified anymore, and that the 
plot in that regime has no other physical significance than 
to provide an important cross check for the DMRG. 

The exact result at (d/ABcs)7 = A describes the two 
grains as a single superconductor with half the level spac- 



FIG. 7: Same calculation as in Fig. but at larger inter- 
grain coupling. Now, the DMRG has a somewhat larger range 
of validity. For small grains, the DMRG curve lies higher 
(i.e. the tight-binding approach is not as good as the DMRG 
anymore) . 



ing d 2 = d/2 and with the interaction Hamiltonian 



H 2 = -X 2 d 2 Y, h \ h r ( 24 ) 

i£R,L, j£R,L 



and with X 2 = 2A. In the large-coupling regime, the 
Josephson energy Ej = £ con d,2 - 2i? cond! i is entirely 
dominated by the condensation energy -E C ond,2 of the 
large superconductor described by Eq. ipijl. which is 
much larger than the condensation energies 2£ , con( j.i of 
the isolated grains (i.e. for 7 = 0). In the BCS limit, 
Ej s» E cond2 = WDcbyc«sinh 2 (l/A 2 ). In particular, Ej 
is seen to be an extensive quantity, i.e. -Ej/wDebyc °c n 
for the particular choice 7 = (Abcs/gOA, for which the 
two superconductors are described as one. In this case 
the inter-grain coupling acts like a bulk term (and no 
longer as a surface effect), which is manifest in the way 
that 7 scales with the system size: 7 scales no longer as 
a constant, but with the volume of the system. 

As is evident from Fig. |S1 the tight-binding method, 
which is only applicable at very small values of 7, ceases 
to be valid long before the point (d/ABCs)7 = A is 
reached. The inset of Fig.|Hlshows an enlargement of the 
main figure for small 7. It is seen that for E® <§; Abcs> 
the results from the tight-binding method and from the 
DMRG agree, as expected. 
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tor and Metallic Clusters". US and DG acknowledge support 
through a Gerhard Hess prize of the Deutsche Forschungsge- 
meinschaft. We thank G. Falci, Y. Imry, S. Kleff, C. Kollath, 
M. Schechter, and G. Sierra for discussions. 
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Tight binding vs. strong coupling approach. X= 0.3, n = 100 
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FIG. 8: The Josephson energy is plotted as a function of 
the inter-grain coupling 7 in grains that are small enough so 
that the DMRG approach works for all values of 7 down to 
zero. Once the coupling is too large (7 ~ 0.06), the tight- 
binding model fails as asserted in Fig. |K| The inset shows an 
enlargement for small values of 7, and illustrates the condition 
Ej <C A sp ~ Abcs for the tight-binding model to be valid, 
which was motivated at the beginning of section ITT1 



APPENDIX A: THE DMRG ALGORITHM IN 
ENERGY SPACE 



In this appendix, some technical aspects of the DMRG 
procedure for approximating the ground state \tp) are ex- 
plained. There are excellent pedago gical reviews of the 
DMRG algorithm to be found [1J,|2(||, as well as a descrip- 
tion of its application to superconducting grains |rl Il7|. 
so we only highlight the key concepts of the DMRG al- 
gorithm for reducing the Hilbert space. Then, the full 
DMRG algorithm as applied in energy space is sketched. 
Finally, a few peculiarities are mentioned that are of rel- 
evance when the algorithm is applied to the problem of 
two coupled superconductors. 

First, we will give an account of the procedure that 
projects out a reduced number of basis states. The 
Hilbert space is divided into two blocks A of states be- 
low and B of states above the Fermi surface, as de- 
picted in Fig. each being represented by the respec- 
tive basis state \i)a and \j)b- A general many-body 
state is expressed as J2ij i'n \^)a ® \j)b- The goal is 
to find a reduced number m of most relevant states 
|u Q ),4 and \up)B, in the sense that they allow for the 
best approximation of the state \tp), such that the norm 

tp) — Ylij rfafiluajA <S> \up)B is minimized, when vari- 
ation over both ij) a p and the states jita)^ \up)B are 
allowed, but only m states per block are to be kept. 
It turns out that the states with this property are pre- 
cisely those eigenstates of the reduced density matrix of 
the respective block (A or B) that correspond to the m 
largest eigenvalues [2fjj. Of course, the larger m is, the 
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FIG. 9: Sketch of the procedure for projecting out the rele- 
vant states in the case of the two-grain DMRG. The shading 
indicates the part of the Hilbert space where only a limited 
number of states are kept. First, a new level is added on grain 
1 (left part of figure). Then, the m most relevant states are 
projected out and kept (right part). Then, a new level on 
grain 2 is added (not shown). 



more accurate the algorithm becomes, until convergence 
is achieved. Typical values for m are m ~ 100 — 400. 

The prescription for the DMRG algorithm is the fol- 
lowing: (i) Start with only a few (2 or 3, say) energy 
levels, few enough that the exact basis of the many-body 
system can be kept explicitly, (ii) Add an additional en- 
ergy level to block A and B, as depicted in Fig.|5]for the 
case of the two-grain DMRG. Construct a basis lua}^ 
for block A, using the basis states from the previous step 
and the exact basis of the newly added energy level. Do 
the same with block B. (iii) Calculate the target state 
ki>), in our case the ground state of the BCS Hamilto- 
nian, within the present Hilbert space, (iv) Calculate 
the reduced density matrix of for block A and B, say 
Pa and pb, by tracing out the full density matrix |V , )(V'| 
over the respective other block. Find the m eigenvec- 
tors |m q )a, \u(3)b, ce,(3 = l.,m, corresponding to the m 
largest eigenvalues of pa and ps- Those are the states 
to be kept as basis states, (v) Transform all operators 
to the new basis. If the blocks A and B are related by a 
symmetry, it may be sufficient to calculate only one set 
of states \u a ). Continue with step (ii) and iterate, until 
the final number of energy levels is reached. 

In step (iii), the ground state \4>) is found using the 
Lanczos procedure, which is very efficient due to the 
sparse nature of the Hamiltonian, but which requires 
many multiplications of a state with the Hamiltonian. 
Since the Hamiltonian is a sparse but extremely large 
matrix (of order m 2 x m 2 ), it is essential not to store 
it as a whole, but to reconstruct it from simple op- 
erators acting only on the blocks A and B when the 
multiplication is performed. For this to be numerically 
possible, it is necessary that the interactions between 
the blocks factorize to a large degree, such that they 
can be expressed as a sum of only a few terms. In 
real-space DMRG, this is always the case as long as 
the interactions are more or less local, but the long- 
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range interactions in energy space do not always factor- 
ize. Luckily, the reduced BCS interaction does factorize 
nicely: H BC s = ~K b a b a + b \ b b + A ^ B), where 
b a,b — J2ieA b c iT c il- A similar factorization is possible 
for the inter-grain coupling IL'l'li in the two-grain DMRG, 
but not for (|T3|) . 

It is also essential for numerical efficiency to make use 
of conserved quantum numbers. In our case, due to par- 
ticle number conservation, it is not necessary to keep all 
the m 2 states |u q )a® \ufj) s as a basis. In our algorithm, 
we keep track of the number l a , lp of particle or hole ex- 
citations associated with each basis vector |m q )a, b, 
respectively. Then, only the states |m q )a <8> | M /3)s have 
to be kept for which 

la - If) = hot, (Al) 

where /tot is the deviation of the total electron pair num- 
ber from half filling. 

In the tight-binding calculation, taking the matrix cl- 
ement (n \bAn + 1) involves approximating two states si- 
multaneously, namely the ground states |n) and n + 1) 
that correspond to the respective number of electron 
pairs n and n + 1. This is simply done by calculating 
both states in step (iii), and by taking the reduced den- 
sity matrix of the mixed state with equal weight in step 
(iv). 

In the two-grain DMRG, the calculations are per- 
formed in the regime that two states, \v) and \v + 1), 
as defined in Eq. (J3J), are degenerate. This is done by 
setting the offset between the energy levels on the left 
and the right grain to zero, and by including one more 
electron pair than there would be at half filling, which 
amounts to setting Ztot = 1 in Eq. I)A1J| . This extra pair 
can, then, be on the left or the right grain at equal energy 



cost. 

One complication arises away from half filling (i.e. 
when Z to t ^ 0): When, in step (iv), the reduced ba- 
sis of one block (block A, say) is calculated by tracing 
over the states in the other block, the part of the trace 
relevant for the states with quantum number l a is, due 
to Eq. IjAljl . performed over states which carry a dif- 
ferent quantum number I p. The dimensionality of the 
two subspaces H(l a ) and T~i(lp) spanned by the part of 
the reduced density matrix with the respective quantum 
numbers might be quite different. However, the rank of 
the reduced density matrix used in step (iv) is limited 
by the dimension of the space over which the trace is 
performed, and therefore, the DMRG only works well as 
long as the dimension oiTi{lfj) is larger than the number 
of states with quantum number l a to be kept. This is 
not guaranteed away from half filling, i.e. when l a ^ lp. 

The problem is solved by mixing a small part (20%) 
into the reduced density matrix that corresponds to the 
ground state at half filling (Z to t = 0). This state will have 
a similar information content as the target state away 
from half filling, as far as the relevant basis vectors are 
concerned, and adds enough to the rank of the reduced 
density matrix for the DMRG to work well. 

In the two-grain DMRG, the energy levels are added 
one by one as depicted in Fig. |5J First levels on grain 1, 
and only afterwards levels on grain 2 are added. They 
are added one by one in order to keep the Hilbert space 
as small as possible. It is also possible and, in fact, would 
be more symmetric, to add both levels at once, but only 
at the cost of having the Hilbert space larger by a factor 
of 4. As it turns out, it is numerically more efficient 
(yielding higher accuracy at the same computation time) 
to add the levels one by one. 
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